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ABSTRACT 

We use a Monte Carlo code to calculate the geodesic orbits of test particles around Kerr black holes, 
generating a distribution function of both bound and unbound populations of dark matter particles. 
From this distribution function, we calculate annihilation rates and observable gamma-ray spectra for 
a few simple dark matter models. The features of these spectra are sensitive to the black hole spin, 
observer inclination, and detailed properties of the dark matter annihilation cross section and density 
profile. Confirming earlier analytic work, we find that for rapidly spinning black holes, the collisional 
Penrose process can reach efficiencies exceeding 600%, leading to a high-energy tail in the annihilation 
spectrum. The high particle density and large proper volume of the region immediately surrounding 
the horizon ensures that the observed flux from these extreme events is non-negligible. 

Keywords: black hole physics - accretion disks - X-rays:binaries 


1. INTRODUCTION 

Prompted by the recent paper bv lBanados et all ([2009( 1 
[BSW], there has been a great deal of interest in the 
potential of Kerr black holes to accelerate particles to 
ultra-relativistic energies and thus to probe a regime of 
physics otherwise inaccessible. The vast majority of this 
work has been analytic and thus largely limited to the 
most simple photon and particle trajectories in the equa¬ 
torial plane. Here we present a more numerical approach 
that focuses on calculating the fully relativistic distribu¬ 
tion function of massive test particles around a spinning 
black hole. With this distribution function and a sim¬ 
ple model for the dark matter annihilation mechanism, 
we can then calculate the annihilation rate and observed 
spectrum as a function of black hole spin and observer 
inclination. 

It has been noted repeatedly in recent works that the 
net energy gained through the Penrose process is quite 
modest, as is the fraction of collision products that might 
escape, and thus the astrophysical imp ortance of the 
BSW effect is quest i onable ((Jacobson fc Sotirioul 120101: 
Banados et all 1201 It lHarada et all l2012t iBeieer et al.l 
20121: IMc Williams! 12013l) . We argue here that two pri¬ 
mary factors (to our knowledge largely neglected in pre¬ 
vious work) could greatly enhance the astrophysical rele¬ 
vance and observability of this annihilation. The first 
is an energy-dependent cross section for dark matter 
(DM) annihilation. This could take many f orms, the 
simp l est of which are p-w a ye annihilation dBertone et al.1 
120051 : iChen fc Zhoull2013t iFerrer fc Hunteidl2013l) . where 
the cross section scales like the relative velocity, or a 
threshold energy, above which the cross section increases 
greatly. This latter assumption is a natural choice for 
a model that includes multiple DM species, with the 
more massive particles intermediate products in the an- 
nihilation process towards gamma rays [see, e.g., IZurekl 
(j2014T) ]. Because gravity is the only known force capa¬ 
ble of accelerating dark matter particles to high energies, 
it is possible that new annihilation channels could occur 
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around black holes that are completely inaccessible ev¬ 
erywhere else in the universe. 

The other effect considered here is the relativistic en¬ 
hancement of the density close to the black hole. This 
is due to the time dilation of observers near the hori¬ 
zon. In a steady-state system, one can think of dropping 
particles into the black hole from infinity at a constant 
rate Poo as measured by coordinate time t. To an ob¬ 
server near the black hole measuring proper time r, an 
enhanced rate T 00 (dt/dr) is seen, with dt/dr > 1. For 
annihilation rates that scale like the density squared, the 
local annihilation rate will be enhanced by (dt/dr) 2 . Of 
course, the products will get redshi fted o n their way back 
out to an observer at infinity (jMcWilliamsll2013lh but we 
are still left with a net enhancement of dt/dr. 

Even without this relativistic enhancement, numer¬ 
ous models also predict an astrophysical enhancement of 
the dark matter density in the galactic nucleus. Adi¬ 
abatic growth of the central black hole will capture 
a large number of particles onto tightly bound orbits, 
growing a steep density spike as the black hole grows 
([Gondolo fc Silkl 119991 : iSadeghian et al.l l201.H l . Gravita¬ 
tional scattering off the dense nucle ar star c lu ster will 
also lead to a dark matter spike ([Gnedin fc Primackl 
l2004f) . similar to the classica l two-body scattering re¬ 
sult of iBahca ll fc Wolj (I1976I ). A t the same time, self- 
annihilation dGondolo fc Silk 1999) a nd elastic s cattering 
((Shapiro fc Paschalidisll2014HFielcls et al.1 12014( 1 will act 
to flatten out this spike into a shallow core more similar 
to the unbound population. 

Because our approach to this problem is predominantly 
numerical, we can easily include treatment of a range 
of black hole spins, particle distributions, and cross sec¬ 
tions, and not limit ourselves to special cases with an¬ 
alytic solutions. Therefore, we can calculate how often 
those extreme cases are likely to occur in a real astrophys¬ 
ical setting [for a notable exception to the analytic ap¬ 
proaches of ear lier work, see th e exha ustive Monte Carlo 
calculations of iWilliama (I1995L 12004( 1 that explored the 
limits of the Penrose process in the context of Comp¬ 
ton scattering and pair production in accretion disks and 
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jets]. Of particular interest has been the following ques¬ 
tion: for two particles each of mass m x falling from rest 
at infinity and colliding near the black hole, what is the 
maximum achievable energy for an escaping photon? We 
find that this limit exceeds 12 ?tj x for an extremal black 
hole with a/M = 1, signific antl y hig her than previously 
published values of 2.6 m x dBeiger et al.l 120121 ). We ex¬ 
plain the underl ying reason forthis discrepancy in a com¬ 
panion paper (lSchnittmanll2014t ). 


2. POPULATING THE DISTRIBUTION FUNCTION 
2 . 1 . Initial conditions 


The primary goal of this paper is to calculate the 8 - 
dimensional phase-space distribution function d/(x, p) of 
DM particles around a Kerr black hole. Two of these 
dimensions are immediately removed due to the assump¬ 
tion of a steady-state solution and stationarity of the 
metric, and the mass-shell constraint of the particle mo¬ 
mentum, leaving us with df(r , 0 , </>,p r ,p 0 ,p 0 ). This func¬ 
tion is further reduced to five dimensions because ax- 
isymmetry removes the dependence on </>. 

To calculate the distribution function, we first distin¬ 
guish between two basic populations: the particles grav¬ 
itationally bound and unbound to the black hole. The 
properties of the bound populations are more sensitive to 
underlying astrophysical assumptions, and will be dis¬ 
cussed below in Section The unbound population 
is more straightforward: we simply assume an isotropic, 
thermal distribution of velocities at a large distance from 
the black hole. Here, “large distance” is taken to be the 
influence radius 7'i n fl of a supermassive black hole with 
mass M, and the DM velocity dispersion is set equal to 
the stellar velocity dispersion ao of the bulge (thus the 
“unbound” population considered in this paper is still 
gravitationally bound to the galaxy, just not_the black 
hole). From the “M-sigma” relation (iFerrarese fe Merrittl 
120001 ) we take 

M « 2 x 1O 7 M 0 ( -(1) 

° V 100 km/s ) K ’ 


and 

GM o ( a 0 V 

rinfl = ~a§~ ~ 8 PC (^Yookin/s ) 


( 2 ) 


In units of gravitational radii r g = GM /c 2 , the influence 

radius is typically quite large: n n fl ~ 10 7 M 7 X ^ 2 r gi where 
M r = (M/1O 7 M 0 ). 

Given this outer boundary condition, we shoot test 
particles towards the black hole with initial velocities 
drawn from an isotropic thermal distribution with char¬ 
acteristic velocity <7o- As we are only interested in 
the distribution function relatively close to the black 
hole, we can ignore any particle with impact param¬ 
eter greater than ss 1000 r g . For those particles that 
we do follow, we calculate their geodesic trajectories 
with the Hamiltonian appr oach described in detail in 
ISchnittman fc Krolikl ( 2013 ) and used in the radiation 
transport code Pandurata. A schematic of this pro¬ 
cedure is shown in Figure |T] As the particle moves 
around the black hole and passes through different fi¬ 
nite volume elements, the discretized distribution func¬ 
tion df(ri,9j, p) is updated with appropriate weights. 


The great advantage of this Hamiltonian approach is 
that the integration variable i s the coordinate time t i n 
Boyer-Lindquist coordinates (iBover fc Lindquist] Il967lf . 
Because this is the time measured by an observer at infin¬ 
ity, it determines the rate at which particles are injected 
into the system in the steady-state limit. Then the distri¬ 
bution function can be populated numerically by assign¬ 
ing a weight to each bin in phase space through which 
the test particle passes, with the weight proportional to 
the amount of time t spent in that volume. The process 
is repeated for many Monte Carlo test particles until the 
5-dimensional distribution function is completely popu¬ 
lated. 


2 . 2 . Geodesics and Tetrads 

Following iSchnittman fc Krolikl (120131 ), we define local 
orthonormal observer frames, or tetrads, at each point in 
the computational volume. Depending on the population 
in question (i.e., bound vs. unbound), it is convenient to 
use eit her the z ero- angular-momentum observer (ZAMO; 
iBardeen et al.l 1 1972 )) or the “free-falling from infinity 
observer” (FFIO) tetr ads. In all cases we use Boyer- 
Lindquist coordinates (IBover fc Lindauisti 119671) . where 
the metric can be written 


9^ = 
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This allows for a relatively simple form for the inverse 
metric: 


<r = 


(-l/a 2 
0 
0 

\—uj/a 2 


0 0 

A/p 2 0 

0 1/p 2 

0 0 


— ui/a 2 

0 

0 

1/w 2 — u} 2 /a 2 


with the following definitions: 
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(5b) 

(5c) 
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Unless explicitly included, we adopt units with G = c = 
1 , so distances and times are often scaled by the black 
hole mass M. 

The ZAMO tetrad can be constructed by 


e (t) “ ~ e W 


— 




J (0)' 
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Figure 1 . Schematic of our method for populating phase space with geodesic trajectories. The test particles are injected at large radius 
(ro = 10 7 r g ) with thermal velocities with dispersion ctq <C c. Those particles passing within 1000 r g of the black hole contribute to the 
tabulated distribution function in each volume element ( ri,Qj ) through which they pass, with a weight proportional to the amount of 
coordinate time t spent in that zone. 



where we designate tetrad basis vectors by fl indices, 
while coordinate bases have normal indices. 

To construct the FFIO tetrad, the time-like basis vec¬ 
tor is given by the 4-velocity correspond¬ 
ing to a geodesic with ut = — 1, ug = = 0, and from 

normalization constraints, 


u r 





1/2 


(7) 


Then the spatial basis vectors are constructed via 
a standard Gram-Schmidt method and aligned roughly 
parallel to the Boyer-Lindquist coordinate bases. 

Any vector can be represented by its components in 
different tetrads via the relation 


u = e (/j) u M = e^)U ,J ’, (8) 

whereby the components are related by a linear transfor¬ 
mation E~: 

u fl = E^u> i , (9a) 

u^ = [E~X^. (9b) 

These u& are the components that we use for the tabu¬ 
lated distribution function^ Because of the normaliza¬ 
tion constraints, we need only store three components of 
the 4-momentum in each spatial volume element, making 
the total dimensionality of the distribution function five: 
two space and three momentum. 

In Pandurata, the geodesics are integrated with 

a variab l e t i me step 5 th o rder Cash-Karp algorithm 
dSchnittman &; Krolikl 1201311 . This technique very nat¬ 
urally matches small time steps to regions of high cur¬ 
vature and thus areas of high resolution in the spatial 


1 While these contravariant indices technically refer to 4- 
velocities, and not 4-momenta, we use the terms interchangeably 
in the locally flat tetrad basis, where most of our calculations take 
place. 


grid. For each time step, a weight proportional to the 
coordinate time spent on that step is added to the dis¬ 
tribution function for that particular volume of phase 
space. Because the particle typically remains within a 
single volume element for many time steps, we find that 
interpolation errors are small. 

The spatial momentum components 7/3* can be posi¬ 
tive or negative and span many orders of magnitude. To 
adequately resolve the phase space and capture the rel¬ 
ativistic effects immediately outside the black hole hori¬ 
zon, we find that on order ~ 10 3 bins are required in 
each dimension. If the entire phase-space volume were 
occupied, this would correspond to an unfeasible quan¬ 
tity of data. Fortunately, this volume is not evenly filled, 
so such a hypothetical 5-dimensional array is in fact ex¬ 
ceedingly sparse. In practice, we are able to use a dy¬ 
namic memory allocation technique that only stores the 
non-zero elements of the distribution function. Yet even 
so, a well-resolved calculation can easily require multiple 
GB of data for a single distribution function, and to ad¬ 
equately sample this phase space requires on the order 
of ~ 10 9 test particles, with each geodesic sampled over 
thousands of time steps. Fortunately, this is a trivially 
parallelizable problem, so it is relatively simple to achieve 
sufficient resolution in a reasonable amount of time with 
a small computer cluster. 

2.3. Unbound Particles 

As mentioned above, for the unbound population, the 
outer boundary condition for the phase space density at 
ring is relatively well-understood. The velocity distribu¬ 
tion is thermal with characteristic speed cr(0- The spatial 
density of dark matter is measured from galactic rotation 

2 While there could be some small anisotropy in the dark matter 
velocity distribution at r; n g, it is unlikely to be correlated with 
the black hole spin. Thus the predominantly radial velocities of 
incoming particles will be independent of polar angle, and therefore 
for all intents and purposes appear isotropic from the black hole’s 
point of view. Similarly, even if the DM velocity distribution at the 
influence radius is not strictly Maxwellian, this too will have little 


















4 


SCHNITTMAN 


curves at kpc distances from the nucleus, and then must 
be extrapolated in to pc distances with a combination 
of observations and stellar profile modeling. For exam¬ 
ple, in the Milky Way the DM density near the Sun is 
0.3 GeV/crn 3 , and the radial profile can be reasonably 
well-modeled with a simple p ~ R~ l profile, giving a 
density of ~ 10 3 GeV/cm 3 at r; n fl. Inside of r- m a there is 
almost certainly an additional bound c omponent to the 
DM distribution dGondolo fc SiUd[i999 .). so the unbound 
population described here can best be understood as a 
strict lower bound on the phase space density. 

Outside of ~ 100 r g the unbound population can be 
treated as a col lisionl ess ga s of accreting particles, as in 
IZeldovich fc Novikov (1971:). In the Newtonian limit, the 
density and velocity dispersion can be written 



r/M 


n(r) = n 0 



2GM \ 1/2 
F ) 


( 10 ) 


and 


\r) = 



2 GM\ 


( 11 ) 


Figure 2. Spatial density (a) and mean relative momentum 
(7/3) rei (6) of unbound particles as measured in the FFIO frame 
in the equatorial plane of a Kerr black hole with a/M = 1. Th e 
dashed lines are the Newtonian solutions of equations (1 101 HID , 
while the solid curves come from the fully relativistic Monte Carlo 
calculation. 


In Figure [2] we show the spatial density of unbound 
particles as measured by a FFIO around a Kerr black 
hole with spin parameter a/M = 1, as well as the mean 
particle momentum as measured in that frame. We find 
very close agreement to the Newtonian results all the 
way down to r ~ 10 r g . The deviation of the momentum 
from the Newtonian solution is due largely to the special 
relativistic terms proportional to the Lorentz boost 7 . 

The proper density is governed by two competing rel¬ 
ativistic effects. One is time dilation and the other is 
spatial curvature. Close to the black hole, the parti¬ 
cle’s proper time r slows down relative to the coordinate 
time t measured by an observer at infinity, giving a large 
dt/dr. This has the effect of increasing the number den¬ 
sity because, in a steady state, particles are injected into 
the system at a constant rate—as measured by an ob¬ 
server at infinity. The injection rate measured by an 
observer close to the black hole is higher by a factor of 
dt/dr , leading to her seeing a larger proper density. 


Figure 3. Proper volume measured in the FFIO frame. The 
dashed line is the Newtonian value dV/dr = Anr" 2 , and the solid 
curve measures the FFIO’s proper volume dV/dr. 




r/M 


r/M 


impact on the results presented here, because the initial velocities 
are so small compared to the orbital velocities near the black hole, 
the trajectories are indistinguishable from particles injected from 
infinity with zero velocity. 


In fact, the proper density would be even higher if 
it weren’t for another important relativistic effect: the 
stretching of space around a black hole. Specifically, 
the Boyer-Lindquist radial coordinate element dr cor¬ 
responds to a greater and greater proper distance as 
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the observer approaches the horizon. This naturally 
gives a greater proper volume dV , shown as a solid 
curve in Figure [3] Again, we show the Newtonian value 
dV/dr = 4nr 2 as a dashed curve. Because the parti¬ 
cle interaction rates scale like n 2 v dV, all these effects 
combine to increase the importance of reactions near the 
black hole. 

In FigureUwe plot the momentum distributions of un¬ 
bound dark matter particles, as observed by a FFIO in 
the equatorial plane, at a relatively large distance from 
the black hole: r = 100M. Each 1-dimensional distribu¬ 
tion is calculated by integrating over the other two mo¬ 
mentum dimensions. We also plot the momentum mag¬ 
nitude 7|/3| in panel (a). Because the particles all have 
relatively small velocities at infinity /3o ~ oo/c <C 1, their 
velocities in the weakly relativistic region r g -C r <C tq 
are given by v « y/2GM/r, corresponding to v « 0.14c 
for r = 100 M. 

For the three spatial components of the momentum 
distribution, we see a nearly isotropic velocity distribu¬ 
tion with a few subtle but interesting deviations. First, 
we note how there is a slight deficit of particles with 
positive p r . This is due to capture by the black hole of 
particles coming in from infinity with nearly radial tra¬ 
jectories. By definition, these particles also have small 
values of p e and p depleting the distribution function in 
those dimensions around /3 = 0. While the distribution 
in the 9 dimension is symmetric, note that the depletion 
in the </> distribution is offset to slightly negative values 
of pA This is due to the well-known preferential capture 
by Kerr black holes of retrograde particles with angular 
momenta aligned opposite to the black hole spin. 

In Figure 0 we plot the phase-space distribution for 
the same boundary conditions as in Figure 01 but now 
at r = 2 M. The difference is quite dramatic, but all the 
features are essentially due to the same physical mecha¬ 
nisms. This close to the horizon, there is a very strong 
depletion of outgoing particles with p r > 0, as most par¬ 
ticles are captured by the black hole. The only particles 
that can avoid capture at this radius have prograde tra¬ 
jectories in the equatorial plane. Thus, the distribution 
is now peaked around p e = 0 instead of showing a deficit. 

There is also a strong peak near p® = 1 due to the 
relatively stable, long-lived prograde orbits that circle 
the black hole multiple times before getting captured or 
escaping back out to infinity. In fact, the distribution 
of coordinate momentum is significantly more lopsided 
to p^ > 0, but this is masked in Figure [SH because this 
distribution is measured by an observer with > 0 
herself. The sharp fall-off of the azimuthal distribution 

above p^ ~ 1 is due to the angular momentum barrier of 
the black hole. Particles with higher values of p^ simply 
never reach this small radius. 

To the best of our knowledge, these distribution func¬ 
tions have never been calculated before for a Kerr black 
hole. However, the particle number density can be de¬ 
termined analytically for a non-spinning Schwarzschild 
black hole, in the limit of oo -C c. This allows at 
least one test of our numerical methods, although admit¬ 
tedly not a very strong one, as most of the interesting 
features are related to the far more complicated orbits 


around a sp inning black hole. We follow the approach 
of lBaushevi (12009H . who integrates the distribution func¬ 
tion with fixed energy, carefully setting the angular mo¬ 
mentum integration bounds based on which orbits are 
captured from a given radius. The results are shown in 
Figure [6l with our numerical calculation plotted as a red 
curve and the analytic result in black, showing perfect 
agreement. Note that Baushev’s expression is given for 
a coordinate density rather than a proper density, which 
also explains the sharper peak at small r. 

2.4. Bound Particles 

As mentioned above, the unbound population can be 
thought of as a lower limit on the total DM density. 
There will also likely be a substantial population of par¬ 
ticles that are gravitationally bound to the black hole. 
As described in lGondolo fc Silkl (|1999H . the origin of the 
bound population is the adiabatic growth of the super- 
massive black hole on a timescale much longer than the 
typical orbital time. This physical mechanism can be un¬ 
derstood as follows: as a marginally unbound DM parti¬ 
cle passes within r; n fl, a small amount of baryonic matter 
is accreted into this region, deepening the potential well 
just enough to capture the particle onto a marginally 
bound orbit. Once captured, the particle continues to 
orbit the black hole while conserving its orbital angu¬ 
lar momentum as the black hole continues to gain mass. 
This has the effect of shrinking the radius of the orbit. 

Over time, more particles are captured and subse¬ 
quently migrate close r to the black hole, b uilding up 
a steep density spike dGondolo fe Silkl [199 91. Inside of 
the inner-most stable circular orbit (ISCO), there is a 
sharp falloff in the density sp ike due to plunging trajec¬ 
tories (jSadeghian et al.l [2013 9. Here we do not attempt 
to solve for the slope of the density spike at large radii 
but leave it as a free parameter, and fix the density 
at the influence radius as for the unbo und population: 
nhmirid (r) — n Q (r/r Q )~ a . Following iGondolo fe SilU 
con, we also allow for the possibility of a density up¬ 
per bound n a nnih due to annihilation losses occurring over 
very long timescales. 

To populate the phase-space distribution for the bound 
population, we follow a similar method as described 
above for the unbound particles, but instead of launch¬ 
ing them from large radius with a limited range of im¬ 
pact parameters, now we launch them in situ with a 
isotropic thermal velocity distribution, as measured by 
a local ZAMO. These particles begin much closer to the 
black hole, so the r elativistic Ma xwell-Jiittner velocity 
distribution is used (|Jiittnerlll91 If) , with the characteris¬ 
tic virial temperature 0(r) = 1/2 [1 — czamo(t)], where 
czamo(t) — 1 is the specific gravitational binding energy 
of the ZAMO. 

Because many of the particles launched close to the 
black hole get captured, we first integrate their trajecto¬ 
ries for a few orbital periods to ensure they are in fact on 
stable orbits. Only then do they contribute to the tabu¬ 
lated distribution function. Additionally, a small fraction 
of the test particles from the tail end of the velocity dis¬ 
tribution will in fact be unbound, and these are similarly 
discarded. 

As with the unbound distribution, for each step along 
its trajectory, the test particle contributes to the phase 
space distribution a small weight proportional to the 
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Figure 4. Momentum distribution of unbound particles observed by a FFIO in the equatorial plane at radius r = 100 M. All particles 
have nearly unitary specific energy at infinity, so the average particle speed is very close to y/2GM/r = VO.02 c (panel a). In panels ( b-d) 
we show the distribution of the individual momentum components, which are nearly isotropic this far from the black hole. 



yIPI 



tP‘ 


amount of time spent on that step. Yet now, instead of 
using the coordinate time dt, we use the proper time of 
the ZAMO frame from which the particles are launched, 
including an additional weight to ensure the appropriate 
radial form of the density distribution at larger radii. 

In Figure [7] we plot the radial density distribution and 
mean relative momentum of the bound particles, as mea¬ 
sured in the ZAMO frame, in the equatorial plane around 
a Kerr black hole with spin a/M = 1. The density pro¬ 
file is constructed so that p(r) ~ r~ 2 at large radii. We 
clearly see major differences relative to the unbound pop¬ 
ulation shown in Figure [2] Because of the lack of sta¬ 
ble orbits close to the black hole, the bound population 
declines inside r ss AM, which corresponds roughly to 
the mean radius of the ISCO for randomly inclined or¬ 
bits around a maximally spinning black hole. This effect 
was described in detai l for non-spinning black holes in 
ISadeehian et al.l (I2013f ). For equatorial circular orbits, 
only prograde trajectories are allowed inside of r = 9 M. 
This leads to all particles moving in roughly the same 
direction closer to the black hole, and explains why the 
relative momentum {'y/3} ie i does not increase nearly as 
fast for the bound population as it does for the unbound 
population, which allows plunging retrograde trajecto¬ 
ries, and thus more “head-on” collisions. 

In Figure [8] we show the 2D density profile in the 
x — z plane for both bound and unbound populations, 



# 





for a/M = 0 and a/M = 1. The horizon in Boyer- 
Lindquist coordinates is plotted as a solid black line. For 
comparison purposes, the density scale is normalized to 
the mean value at r = 10M. In reality, the density of 
the bound par ticles could be orders o f magnitude greater 
at these radii dGondolo fe Silkl 117) 99). The most obvious 
difference here is the depletion of bound orbits inside of 
the ISCO, which lies at r = 6 M for non-spinning black 
holes. For spinning black holes, the radius of the ISCO 
is a function of the particle’s inclination angle, ranging 
from r = 1 M for prograde orbits in the equatorial plane, 
to r = 5.2 M for polar orbits, and r = 9 M for retrograde 
equatorial orbits. 

Inside of the ISCO, there is also the “marginally 
bound” radius, where particles with unity specific energy 
can exist on unstable circular orbits. This radius is also 
a function of inclination angle, and is plotted in Figure [8] 
as dotted curve. Inside of this orbit, no bound particles 
will be found (for improved visibility, we have left this 
region white, not black, as would be required by a strict 
adherence to the color scale). One interesting feature of 
Figure [8] is that the density of the unbound population 
around spinning black holes doesn’t show any obvious 
^-dependence. It appears that the enhanced density due 
to long-lived prograde orbits is almost exactly countered 
by the lack of retrograde orbits at the same latitude. 

In Figure [9] we show the phase space distribution 
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Figure 5. Momentum distribution of unbound particles observed by a FFIO in the equatorial plane at radius r = 2 M. Unlike Figure HI 
here we see a decidedly non-thermal and highly anisotropic distribution. 
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Figure 6. Comparison of our numerical results (red) with the ana¬ 
l ytic e xpression (black) for the particle density derived bv IBaushevI 
120091 ') for a Schwarzschild black hole. The density here is defined 
in the coordinate, not proper, frame, leading to a much steeper 
rise at small r. In Boyer-Lindquist coordinates, the horizon for a 
non-spinning black hole is at r = 2 M. 



r/M 


for each of the momentum components, as measured 
by a ZAMO in the equatorial plane at large radius 
(r = 100M). Compared to the equivalent plot for the 
unbound distribution (Fig. [4]), we see a number of sig¬ 
nificant differences. First, the fact that these particles 



# 





are bound requires that E < 1, and the imposed virial 
energy distribution results in mean velocities that are 
smaller than those of the unbound population by a fac¬ 
tor of ~ y/2. Second, because we require stable, long- 
lived orbits, there is a larger depletion around p 9 = 0 
and p$ = 0, as these trajectories are all captured by the 
black hole and thus do not contribute at all to the dis¬ 
tribution function. Similarly, we see a larger asymmetry 
due to the preferential capture of retrograde orbits with 

p$ < 0 . 

In Figure [TU] we plot the same momentum distribution 
functions, now at r = 2 M. Here the contrast with the 
unbound population (Fig. [5}) is even greater. The only 
stable orbits at this radius are prograde, nearly circular, 
nearly equatorial orbits. This results in a relatively nar¬ 
row distribution clustered around u& = [%/2,0, 0,1] in the 
ZAMO frame. This narrower range in allowed velocities 
will have a profound impact on the shape of the annihi¬ 
lation spectrum, as we will see in the following section. 


3. ANNIHILATION PRODUCTS 

Once we have populated the distribution function, 
we can calculate the annihilation rate given a simple 
particle-physics model for the dark matter cross section. 
Again, it is simplest to work in the local tetrad f rame . 
Including special relativistic corrections (| Weaved 1197611 . 
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Figure 7. Spatial density (a) and mean relative momentum 
(7/3) rei (6) of bound particles in the equatorial plane of a Kerr 
black hole with a/M = 1, as measured in the ZAMO frame. The 
dashed lines are the Newtonian solutions when n(r) ~ r —2 far from 
the black hole. 
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the local reaction rate is given by the following: 

i?(x) = [ e? 3 pi [ d 3 p2 /(x, pi) /(x, p 2 ) TfL (j x (7 re i) u re i, 
J J 7 i 72 

( 12 ) 

where 71 and 72 are the Lorentz factors of two parti¬ 
cles as measured in the tetrad frame, u re i is their relative 
velocity, and a x is the annihilation cross section (poten¬ 
tially a function of the relative velocity), -ff(x) has units 
of [events per unit proper volume per unit proper time], 
so we multiply by dr/dt to get the rate observed by a 
distant observer. 

The distribution function /(x, p) is calculated numeri¬ 
cally using the methods of Section [2] As discussed there, 
the numerical representation of / can have upwards of 
10 s elements, so the direct integration of equation m 
is generally not computationally feasible. Instead, we use 
a Monte Carlo sampling algorithm to pick random mo¬ 
menta for each particle with an appropriate weight based 
on the magnitude of / and the size of the discrete phase 
space volume. 

The spatial integration, however, is carried out di¬ 
rectly, looping over coordinates r and 9. This is shown 
schematically in Figure [11] For each volume element, a 
large number (typically ~ 10 6 ) of pairs of particles are 
sampled, and for each pair, a center-of-mass tetrad is 


created. The total energy in the center-of-mass frame is 
given by 

^com = m%\J 2(1 + U! • U 2 ) , (13) 

where m x is the rest mass of the DM particle, and u = 
p /m x is the particle 4-velocity. 

The 4-velocity of the center-of-mass frame is then given 
by 

Ucom = (ui + U.2)/£'com • (14) 

The center-of-mass tetrad is constructed with = 
Ucom- The spatial basis vectors are totally arbitrary, as 
they are only needed to launch photons with an isotropic 
distribution in the center-of-mass frame. Two photons, 
labeled k 3 and k 4 in Figure [TT] are launched in opposite 
directions, each with energy E com /2 in the center-of-mass 
frame. We then transform back to a coordinate basis for 
the geodesic integration of the photon trajectories to a 
distant observe r. 

As in iSchnittman fc Krolikl (120131) . for the photons 
that reach infinity, Pandurata can generate an image and 
spectrum of the emission region. An example in shown 
in Figure[l2]for the annihilation signal from the unbound 
population around an extremal black hole, limiting the 
emission signal to the region r < 100M. While the flux 
clearly increases towards the center of the image, because 
the density and velocity profiles are relatively shallow 
(see Fig. [2] above), the net flux is actually dominated by 
emission from large radii. These annihilation events are 
not very relativistic, so produce a strong, narrow peak 
in the observed spectrum, centered at the DM rest mass 
energy. 

The annihilation events occurring closer to the hori¬ 
zon sample a much more energetic population of parti¬ 
cles. Restricting ourselves to only those events where the 
center-of-mass energy is greater than 1.5 x the combined 
rest mass of the annihilating particles, we can zoom in 
to the center of Figure ITS1 The result is shown in Figure 
m now focusing on the inner region within r < 6 M. At 
these small radii, the effects of black hole spin become 
much more evident. One such effect is the characteristic 
shape of the Kerr shadow, defin ed by the impact pa - 
rameter of critical photon orbits (Chandrasekhar 1983). 
The observed flux is clearly asymmetric, as the prograde 
photons originating from the left side of the image have 
a much greater chance of escaping the ergosphere and 
reaching a distant observer. 

There is another interesting feature of Figure [13] that 
we believe is novel to this work. Namely, the purple lobes 
emerging from the “mid-latitude” regions near the center 
of the image. These are regions of greater photon flux, 
albeit very highly redsliifted. Recall, this image is cre¬ 
ated by considering only annihilations with moderately 
high center-of-mass energy. Near the equatorial plane, 
extreme frame dragging ensures that the velocity disper¬ 
sion is highly anisotropic, with most of the DM particles 
and their annihilation photons getting swept along on 
prograde, equatorial orbits. Above and below the plane, 
the DM distribution is more isotropic, leading to a more 
isotropic distribution of outgoing photons. Yet if one 
goes two far off the midplane, it becomes more difficult 
for the photons to escape. At the mid-latitudes, there is 
just enough frame dragging for photons to escape, yet not 
so much that they get deflected away from the observer. 
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Figure 8. Spatial density of test particles in the x — z plane, for both bound and unbound populations, for a/M = 0 and a/M = 1. For 
each case, we show the unbound distribution on the left side and the bound distribution on the right side of the plot, and all distribution 
functions are normalized to the mean density at r = 10 M. The horizon is plotted as a solid curve and the radius of the marginally bound 
orbits is shown as a dotted curve. The spin axis of the black hole is in the +z direction. 




The spectrum corresponding to this image is also plot¬ 
ted in Figure [13] Not surprisingly, the red and blue 
wings of the annihilation line shown in Figure QT] come 
from the most r e lativis tic events. As pointed out by 
iPiran fc Shahaml (119771 ). even reactions with very high 
center-of-mass energies will typically lead to photons 
with low energies as measured at infinity, thus explain¬ 
ing the red tail of the annihilation spectrum. The high- 
energy tail above E = 2 m x is due exclusively to Penrose- 
process reactions where one of the annihilation photons 
has negative energy and gets cap tured by the black hole 
(|Penroselll969l : IPiran et~allll975h . 

Earlier analytic work predicted that the maximum 
energy attainable from the collisional Penrose process 
was 2.6m v for particles falling from rest at infinity 
(jHarada et al.ll2012LiBeiger et al.l[2012l ). Because our cal¬ 
culation is fully numerical, it was able to reveal previ¬ 
ously unknown trajectories leading to very high efficien¬ 
cies with E > 10 m x , as seen in Figure [TS] Closer inspec¬ 
tion revealed that these high-energy photons are created 
when an infalling retrograde particle collides with a out¬ 
going prograde particle that has just enough angular mo¬ 
mentum to reflect off the centrifugal barrier, providing 
the necessary energy and momentum for the annihila¬ 
tion photon to e scape the black hole (iSchruttmanl 120141 : 
iBerti et aI1l2014D . 

Due to the strong forward-beaming effects within 
the ergosphere, the escaping photon flux is highly 
anisotropic, with the peak flux and highest-energy pho¬ 
tons emitted in the equatorial plane. Figure QT] shows the 
predicted annihilation spectra for observers at different 
inclination angles for the same DM profile as shown in 
Figure [13] Again, we restrict ourselves to the highest- 
energy reactions with E com > 3 m x . 

It is also instructive to plot the annihilation flux as a 
function of the emission radius. In Figure [15] we show 
both the observed flux (solid curves) and the flux that 
gets captured by the black hole (dashed curves) as a func¬ 
tion of radius, integrated over all observing angles. The 
emission is further subdivided by the center-of-mass en¬ 
ergy of the annihilating particles. Of course, the photons 
emitted closer to the black hole have a greater chance of 


getting captured. For the unbound population, the to¬ 
tal escape fraction ranges from f esc = 93% at r = 10M 
down to / esc (2M) = 14%, and / esc (l.lM) = 0.25%. At 
small radius, these numbers are somewh at smaller than 
those calculated bv lBanados et al.l (1201 HI . who only con¬ 
sidered critical trajectories in the equatorial plane, where 
the escape probability is greatest. Yet at large radius, 
our distribution includes particles with typically greater 
impact parameters, and thus greater chance for escape. 

Another interesting feature of the curves in Figure [15] 
is the very sharp cutoff above a critical radius for each 
energy bin. This is a natural consequence of conservation 
of energy. Because all unbound particles come in from 
rest at infinity with E = m x , the available kinetic energy 
in the center-of-mass frame is simply the gravitational 
potential energy Mm x /r at that radius. For example, 
to reach a center-of-mass energy of 10% above the rest 
mass energy, the particles must fall within r ss 10 M. 
Also note that inside r k, AM, most of the photons are 
captured, while outside of this radius, most escape. This 
is in close agreement with what we found for plunging 
orbits inside of the ISCO of a Schwarzschild accretion 
flow in iSchnittman et al.l (120131) . 

On the other hand, for the bound population of DM 
particles (Fig. 115b ). which by definition are not plunging, 
we find that the photon escape fraction is more than 90% 
at all radii, greatly increasing the relativistic effects ob¬ 
servable from infinity. This is consistent with the classic 
calculation bv iThornel (Il974l) which found that for thin 
accretion disks limited to circular, planar orbits outside 
the ISCO, the fraction of emission ultimately captured by 
the black hole was never more than a few percent, even 
for maximally spinning black holes where the majority 
of the flux emerges from extre mely close to the horizon. 

As we showed in iSchnittmanl (I2014T) . the peak en¬ 
ergy attainable from particles falling in from infinity is 
a strong function of the black hole spin. Now, consid¬ 
ering the full phase-space distribution function of the 
particles, we can see how the shape of the spectrum de¬ 
pends on spin. In Figure [TO] we plot the flux seen by 
an equatorial observer, again limited to the high-energy 
annihilations with E com > 3 m x . For even marginally 
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Figure 9. Momentum distribution of bound particles observed by a ZAMO in the equatorial plane at radius r = 100 M. Unlike the 
unbound distribution in Figure [4] the energy distribution is much broader here, yet with a smaller mean momentum (panel a). In panels 
( b-d ) we show the distribution of the individual momentum components. 
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sub-extremal spins, the peak photon energy falls precipi¬ 
tously. As the spin decreases further, the number of colli¬ 
sions with E com > 3 m x also decreases, thereby reducing 
the total flux observed. Lastly, the decreasing spin also 
increases the critical impact parameter for capturing pro¬ 
grade photons, making it harder for the annihilation flux 
to escape to infinity. 

Recall from Section m above that the density of the 
unbound distribution scales like n ~ r -1 / 2 . From the 
rate calculation in equation (1121) we see that the annihi¬ 
lation rate [events/s/cm 3 ] scales like R(r) ~ r _3//2 . In¬ 
cluding the volume factor dV = 4nr 2 dr we can write 
the differential annihilation rate as dR/dr ~ r 1 / 2 . In 
other words, the unbound contribution to the annihila¬ 
tion signal diverges at large radius. In practice, the outer 
boundary can be set as the black hole’s influence radius, 
typically 10 6 ~ 7 r g . This means that the observed signal 
will essentially be a delta function in energy, with only 
small perturbations from the relativistic contributions at 
small r, and thus measuring spin from annihilation lines 
would be a very challenging prospect indeed. 

Two possible effects provide a way around this prob¬ 
lem, each with its own additional uncertainties. One pos¬ 
sibility is that the annihilation cross section is a strong 
function of energy, increasing sharply above some thresh¬ 
old energy. This is admittedly rather speculative, and 
in conflict with leading DM models of self-annihilation 


dBertone et all 1200511 . On the other hand, we do not 
even know what the dark matter particle is, or if there 
are many DM species making up a rich “dark sector,” 
with all t he bea ut y and complexity of the standard model 
particles (iZurekl[20141 1. One could easily imagine a DM 
analog of pion production via the collision of high-energy 
protons, in which case the only reactions could occur im¬ 
mediately surrounding a black hole, the ultimate gravita¬ 
tional particle accelerator. In this case, by construction 
the annihilation rate is dominated the region immedi¬ 
ately surrounding the black hole. 

Another possibility is that the DM density is domi¬ 
nated by a population of bound particles. As described 
above in section m this population arises through 
the adiabatic growth of the black hole through accre¬ 
tion, capturing marginally unbound particles while also 
making the bo u nd pa r ticles ever more t ightly bound 
iGondolo fc SilU (|1999lh iSadeghian et all (|2013h . This 
process will generally lead to a much steeper density pro¬ 
file, such as the n ~ r -2 distribution we use here. In this 
case, the differential reaction rate scales like dR/dr ~ 
r -5 / 2 so the annihilation spectrum is now dominated by 
the particles at smallest radii. In both cases—energy- 
dependent cross sections and a large boun d po pulation— 
the relativistic effects described in Section liOl Cexpanded 
proper volume and time dilation) push the most impor¬ 
tant interaction region to even smaller radii, and thus 
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Figure 10. Momentum distribution of bound particles observed by a ZAMO in the equatorial plane at radius r = 2 M. Compared to 
Figure l5l here we actually see a more symmetric, thermal distribution making up a thick torus of stable, roughly circular orbits near the 
equatorial plane. 
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Figure 11. For a given phase-space distribution /(x, p), the an¬ 
nihilation rate is calculated in each discrete volume element around 
the black hole. Every annihilation event samples the distribution 
function to get the momenta for the two dark matter particles pi 
and p 2 and produces two photons k 3 and k 4 with isotropic dis¬ 
tribution in the center-of-mass frame. The product photons then 
propagate along geodesics until they reach a distant observer or 
get captured by the black hole. 


the annihilation spectra are even more sensitive to the 
black hole spin. 

In Figure [171 we show the annihilation spectra for both 
the bound and unbound populations for a variety of 


spins, now including emission out to r = 1000 M. The 
relative amplitudes are somewhat arbitrary, because we 
don’t know what the relative densities of the two pop¬ 
ulations might be (see discussion below in Sec. SI), but 
it is almost certain that the bound population should 
dominate, possibly ev en by many orders of magnitude 
dGondolo fc Silklfl999tl . At the same time, the unbound 
signal will be even narrower and have a greater ampli¬ 
tude peak than shown here, as it is dominated by low- 
velocity particles at large radius. So while their over¬ 
all amplitudes are uncertain, the detailed shapes of the 
spectra away from the central peak are relatively robust, 
depending only on the properties of geodesic orbits near 
the black hole. 

In this broad part of the spectrum, the bound and 
unbound signals show very different behavior. For non¬ 
spinning black holes, no particle can remain on a bound 
orbit inside of r = 4 M (see Fig. El), so there are no an¬ 
nihilation photons coming from just outside the horizon, 
and these are the photons that produce the most strongly 
redshifted tail of the spectrum. As the spin increases and 
the ISCO moves to smaller and smaller radii, the line be¬ 
comes steadily broader. On the other hand, the unbound 
particles are found all the way down to the horizon, where 
they can annihilate to highly redshifted photons regard¬ 
less of the black hole spin. 

Comparing Figures 151 and flOl we see that the unbound 
particles probe a much greater volume of momentum 
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Figure 12. Simulated image and spectrum of the annihilation 
signal from unbound dark matter out to a radius r = 100M around 
a Kerr black hole. The observer is located in the equatorial plane. 
While the brightness peaks towards the black hole, the total flux 
is dominated by annihilations at large radii. The central shadow 
is clearly seen, blocking emission coming from the far side of the 
black hole. The photon energy E is scaled to the dark matter rest 
mass m x . 
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Figure 13. Simulated image and of the annihilation signal around 
an extremal Kerr black hole, now considering only annihilations 
with Bcom > 3 m x . The observer is located in the equatorial 
plane with the spin axis pointing up. While the image appears 
off-centered, it is actually aligned with the coordinate origin. The 
photon energy E is scaled to the dark matter rest mass m x . 
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space at small radii. This in turn leads to a greater 
chance of producing the extreme Penrose particles that 
characterize the blue tail of the spectrum. Because all 
the bound particles are essentially on the same prograde, 
equatorial orbits, it is much more difficult to achieve 
annihilations with large center-of-mass energies, so the 
high-energy cutoff in the spectrum is much closer to the 
classical result for a singl e particle d ecaying into two pho¬ 
tons in the ergosnhere (|Waldlll974l h In short, for bound 
particles the red tail of the spectrum is a better probe of 
black hole spin, while for the unbound population, the 
blue tail is the more sensitive feature. But in both cases, 
higher spin leads to a broader annihilation line. 

4. OBSERVABILITY 

In addition to the dependence on the dark matter 
density profile, the amplitude of the annihilation spec¬ 
trum will also depend on the unknown dark matter mass 
and annihilation cross section. At this point, it is only 
possible to use existing observations to set upper lim¬ 
its on these unknown parameters. One major obstacle 
that has plagued nearly all observational efforts to detect 
dark matter annihilation is the existence of more conven¬ 
tional astrophysical objects such as active galactic nuclei 


Figure 14. Observed annihilation spectrum for the unbound DM 
population, as a function of observer inclination angle, considering 
only annihilations with E com > 3 m x . The black hole spin is a/M = 
1. 
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(AGN), pulsars, and supernova remnants, all of which are 
powerful sources of high energy gamma rays. One solu¬ 
tion to this problem is to focus on nearby dwarf galaxies, 
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Figure 15. Flux reaching infinity (solid curves) and getting cap¬ 
tured by the black hole (dashed curves), as a function of the center- 
of-mass energy and radius of annihilation, for both bound and un¬ 
bound populations. The black hole spin is maximal. Note that the 
scale on the y-axis is arbitrary, and depends strongly on the anni¬ 
hilation cross sections and peak density. The radial flux profile, on 
the other hand, is a robust result for these populations. 
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Figure 16. Observed spectrum as a function of black hole spin, for 
an observer at inclination i = 90°, considering only annihilations 
with .Ecom > 3 m x . 
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which are thought to have a high DM fraction and are 
not typically conta minated by AGN activi ty or signifi¬ 
cant star formation (lAckermann et alJl201 ill (note, how- 


Figure 17. Comparison of annihilation spectra from bound and 
unbound populations, including all emission out to r = 1000M. 
The peak of the unbound signal will actually be even narrower, as 
it is dominated by annihilations at large radii with small relative 
velocities. 
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ever, the recent work bv iGonzalez-Morales et all ( 2014 ). 
which focuses on the contribution of black holes in dwarf 
galaxies). 

Yet for our purposes, it turns out that the strongest 
upper limits actually come from the most massive galax¬ 
ies with the most massive central black holes. Massive 
elliptical galaxies have the added advantage of being rel¬ 
atively quiesc ent both in nucle a r acti vity and star for¬ 
mation [e.g., iSchawinski et al.l (120071) ]. As mentioned 
above, the annihilation signal from the unbound pop¬ 
ulation will be dominated by flux at large radius. It is 
difficult enough to spatially resolve even nearby black 
holes’ influence radii with HST, much less ganrma-ray 
telescopes, so any potential annihilation signal will tell 
us little about the black hole itself. 

Prospects for detection of an unambiguous black hole 
signature improve if we consider annihilation models 
that include an energy dependence to the dark matter 
cross section. For example, p-wave annihilation mecha¬ 
nisms will have cross sections proportional to the rela- 
tive velocity between the two annihilating pa rticles [see 
iChen fe Zhoul (120131 ) ; iFerrer fe Hunted (|2Q13fl and refer¬ 
ences therein]. Unfortunately, from equation (lllll we see 
that this would only lead to an additional factor of r" 1 / 2 
in the integrand of equation m, which would still be 
dominated by the contributions from large r. 
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Figure 18. Comparison of annihilation spectra from unbound 
populations, for two simple models of the dark matter cross sec¬ 
tion. All spectra are normalized to their peak intensity. For this 
comparison, all emission within r = 10 4 M is included. 
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This effect is shown in Figure |TH1 which plots the pre¬ 
dicted spectra for two annihilation models: v x (v) = 
const (black curves) and u x (v) oc v (red curves). The 
black hole spins considered are a/M = 0 (dashed curves) 
and a/M = 1 (solid curves), and in all cases only the 
unbound population is included. Integrating out to 
r = 10 4 M, we see only a slight difference in the shape 
of the spectrum, with the a x (v) oc v model leading to a 
slightly broader peak (all curves are normalized to give 
a peak amplitude of unity). 

Another possible annihilation model is based on a res¬ 
onant reaction at so me en e rgy ab ove the DM rest mass, 
as suggested in iBaushevI (2009). If the cross section 
increases sharply around a given center-of-mass energy, 
this would have the effect of focusing in on a relatively 
narrow volume of physical space around the black hole, 
as in Figure [15] 

Alternatively, the cross section could abruptly increase 
above a certain threshold energy, if new particles in the 
dark sector become energetically allowed, analogous to 
pion production via proton scattering. In either the 
resonant or threshold models for the annihilation cross 
section, one might imagine a pair of heavier, intermedi¬ 
ate dark particles getting created and then annihilating 
to two photons as in the direct annihilation model. If, 
for example, the mass of these intermediate particles is 
1.5ro x , then the observed spectrum would look like those 
plotted in Figures fl3l and [TGl With a significant increase 
in the cross section above such an energy threshold, these 
relativistically-broadened spectra could in fact dominate 
over the narrow line component produced by the rest of 
the galaxy. 

A less exotic option would be the simple density en¬ 
hancement due to the bound population. If this is suf¬ 
ficiently large, it would easily dominate over the rest of 
the galaxy and also produce a characteristically broad¬ 
ened line sensitive to both black hole spin magnitude and 
orientation relative to the observer. Somewhat ironically, 
one of the things that could ultimately limit the strength 
of the annihilation signal from bound dark matter is an¬ 
nihilation itself. If the adiabatic black hole growth oc¬ 
curred at high redshift, then in the subsequent ~ 10 10 


years, the bound population will get depleted via self- 
annihilation at an acce l erated pace due to its high den sity 
(|Gondolo fc SilklH999t iGonzalez-Morales et alJl2014lf . 

On the other hand, if the black hole grows through 
mergers, or experiences even a single merger since the 
last extended accretion episode, it is quite likely that the 
bound dark matter population could get completely dis¬ 
rupted. The details of such an event are beyond the scope 
of this paper, but could be modeled by following test par¬ 
ticles bound to each black hole through the_merger, via 
post-Newtonian c alculations dSch nittmanl 1201011 or nu¬ 
merical relativity (ivan Meter et al.ll2l)lll i. 

The observational challenge is readily apparent: the 
black holes with the largest bound populations will tend 
to be in gas-rich galaxies with a lot of accretion and 
high-energy nuclear activity that could overwhelm the 
DM annihilation signal. The more massive black holes, 
residing in gas-poor quiescent galaxies, are also more 
likely to have lost their cloud of bound dark matter 
through a history of mergers. Even in the event that 
a gas-rich spiral galaxy hosts a quiescent nucleus, the 
black holes in those ga laxies tend to have lower masses 
dKormendv fc H (ill 201. 'll) . 

While the relation between black hole mass and dark 
matter density is quite complicated for the bound popu¬ 
lation, it is relatively straightforward to calculate for the 
unbound population, which we can take as a lower bound 
on the DM density. Recall the influence radius ?'infl is the 
distance within which the gravitational potential is dom¬ 
inated by the black hole, as opposed to the nuclear star 
cluster or dark matter halo. From equation f2]) we see 
that the influence volume scales like r? nfl ~ M 3 / 2 , while 
the total mass enclosed is—by definition—of the order of 
M. If the dark matter and baryonic matter have similar 
profiles (by no means a certainty!), then more massive 
black holes should have lower surrounding DM density, 
with rii n fl ~ M -1 / 2 . 

Because the unbound DM density falls off more rapidly 
outside the central core, the annihilation flux F un bound 
will be dominated by the contribution from around r; n fl, 
so we can estimate 

^unbound ~ -Jyi nf nfl o- x Ui nfl rf nfl ~ M 3/4 D~ 2 , (15) 

with D the distance to the black hole and the mean ve¬ 
locity at the influence radius = ctq. 

If we consider a threshold energy annihilation model 
where all the flux comes from inside a critical radius 
r cr it ~ few x r g , then the density scales like n cr i t ~ 
riinfl(nnfl/fcrit) 1 / 2 ~ M -3 / 4 while the relative velocity 
scales like v cr jt ~ o’o( ? ’infl/r'crit) 1 / 2 ~ M°. The net flux 
then scales like 

Funbound ~ ~Jj2 G:ritb"xI'V'nt.l" cr it ^ Af ^ 1) . (16) 

In both cases, it appears that the brightest sources will 
be the closest, as opposed to the most massive. 

Now consider the case where the annihilation signal is 
dominated by the bound contribution, the bound den- 
sity is in turn limited by a self-annihilation ceiling as in 
iGondolo fe SilU (|1999f l. and there is a threshold energy 
above which the cross section greatly increases. In this 
case, the flux is simply proportional to the total volume 
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within the critical radius, so Abound ~ M 3 D~ 2 . With 
this scaling, the greatest flux will actually come from 
more distant, more massive black holes. For example, 
NGC 1277, with a mass of 1 .7 x 1 0 lo Mm a nd at a dis¬ 
tance of 20 Mpc (Ivan den Bosch et al.ll2012l ). could give 
an observed flux over a thousand times greater than our 
own Sgr A*! 

Recent works by [Fields et al.1 (I2014T) and 
iGonzalez-Morales et al.1 (l2014f) have argued that 
current Fermi limits of gamma-ray flux from Sgr A* 
and nearby dwarf galaxies with massive black holes 
already place the strongest limits on annihilation from 
DM density spikes. Based on the arguments above, we 
believe that even stronger limits should come from more 
distant, massive galaxies. The other important advance 
presented in the present work is that, for either the 
energy-dependent cross sections, or the steep density 
spikes, the annihilation signal will be dominated by 
the region closest to the black hole, and thus a fully 
numerical, relativistic rate calculation is absolutely 
essential. 

Lastly, we should mention that gamma-rays, while the 
primary observable feature explored in this work, are not 
the only promising annihilation product. High-energy 
neutrinos could also be produced in some annihilation 
channels, particularly those with energ y-dependent cros s 
sections like p-wave annihilation dBertone et al.1 [20051. 
While neutrinos obviously present many new detection 
challenges, the successful commissioning of new astro¬ 
nomical observatories like IceC ube make this approach 
an exciting prospect (lAartsen et al.ll20l3l ). Furthermore, 
the non-DM backgrounds may contribute significantly 
less confusion in the neutrino sky. 

5. DISCUSSION 

As apparent in the previous section, there are still far 
too many unknown model parameters to allow for quan¬ 
titative predictions of the annihilation flux from dark 
matter around black holes. iSadeghian et al.1 ( 2013 ) put 
it best: “There are uncertainties in all aspects of these 
models. However one thing is certain: if the central 
black hole Sgr A* is a rotating Kerr black hole and if 
general relativity is correct, its external geometry is pre¬ 
cisely known. It therefore makes sense to make use of 
this certainty as much as possible.” We have attempted 
to follow their advice to the best of our ability. 

Thus, in order of decreasing confidence, the results in 
this paper can be summarized by the following: 

• For a given DM density rii n fl and velocity dispersion 
(To at the black hole’s influence radius, the fully 
relativistic, 5-dimensional phase-space distribution 
has been calculated exactly for any black hole spin 
parameter, covering the region from rj n fl all the way 
down to the horizon. 

• Given this distribution function and a model for 
dark matter annihilation, the observed gamma-ray 
spectrum can be calculated by following photons 
from their creation until they are either captured 
by the black hole or reach the observer. Two im¬ 
portant relativistic effects serve to increase the an¬ 
nihilation rate as compared to a purely Newtonian 
treatment: time dilation near the black hole effec¬ 
tively raises the density of the unbound population 


in a steady-state distribution being fed from infin¬ 
ity; and transforming from coordinate to proper 
distances greatly increases the interaction volume 
in the region immediately around the black hole 
(see Fig. [3|. 

• Our numerical approach has unveiled previously 
overlooked orbits that can produce annihilation 
photons with extreme energies, far exceeding pre¬ 
vious estimates for the maximu m efficiency of the 
collisional Penrose process (lSchnittmanl[20l4) . The 
peak energy attainable for escaping photons is a 
strong function of the black hole spin. 

• The population of bound dark matter has also 
been calculated numerically, although this depends 
on two additional physical assumptions: a local 
isothermal velocity distribution with a virial-like 
temperature; and an overall radial power-law for 
the density, as found inlGondolo fc~S ilk (|1999|f and 
ISadeghian et al.1 (2013). Including only the long- 
lived stable orbits, we found that the density peaks 
in the equatorial plane somewhat outside of the 
ISCO, forming a thick, co-rotating torus around 
the black hole spin axis. Because the bound pop¬ 
ulation is not plunging towards the horizon, the 
emerging flux has a much greater chance of escap¬ 
ing the black hole. 

• The annihilation spectra from both the bound and 
unbound populations are sensitive to the spin pa¬ 
rameter, but in opposite ways: the unbound spec¬ 
trum varies mostly in the high-energy cutoff, with 
higher spins allowing higher-energy annihilation 
products; the bound population moves closer and 
closer to the horizon with increasing spin, giving 
a stronger red-shifted tail to the annihilation spec¬ 
trum. Both bound and unbound spectra become 
more sensitive to observer inclination with increas¬ 
ing spin, as the spherical symmetry of the system 
is broken. 

• For dark matter particle physics models with an 
energy-dependent cross section (particularly one 
that increases with center-of-mass energy), the an¬ 
nihilation spectrum will be a more sensitive probe 
of the black hole properties. For DM models incor¬ 
porating a rich population of dark sector species, 
black holes may be the most promising way to ac¬ 
celerate these particles and observe their interac¬ 
tions. 

• The shape of the annihilation spectra is relatively 
robust, but the normalization is highly dependent 
on uncertain parameters such as the dark matter 
density profile and cross section. If the unbound 
density profile follows the baryonic matter, with 
the shallow slopes seen in core galaxies, the ob¬ 
served flux should be a relatively weak function of 
black hole mass. If, on the other hand, the anni¬ 
hilation signal is produced by the most relativistic 
population within r cr it ~ few x r g , then the signal 
could scale like M 3 and thus be dominated by the 
most massive black holes in the local Universe. 
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While this paper has treated the bound and un¬ 
bound particles separately, future work will also con- 
sider the self-interaction between these two pop ulations 
dShapiro fc Paschalidisl 120141 : iFields et alJl2014H . which 
may lead to a single, self-consistent steady-state distri¬ 
bution with density slope between —1/2 and —2. Future 
work will also focus on developing a robust framework in 
which we can use existing and future gamma-ray obser¬ 
vations to constrain various parameters of the particle 
physics (e.g., m x , a x (E), and the annihilation mecha¬ 
nism, i.e., line vs continuum) and astrophysical models 
(fiinfl, the bound distribution normalization and slope, 
and the black hole mass, spin, and inclination). While 
initial work will focus on setting upper limits on reac¬ 
tion rates by looking at quiescent galaxies, our ultimate 
ambition is nothing short of an unambiguous detection 
of dark matter annihilation around supermassive black 
holes. 
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